clear;
load('options_BE.mat')
load('Data_FR_revenue.mat', 'Data_FR_revenue')
nsect = length(Data_FR_revenue);
%nsect=1
kjInd_FR_rev = zeros(nsect,1);

for s = 1:nsect
    sigma=Data_FR_revenue(s,3);    
    %sigma=1.0
    k_cool = @(k_j)sd_log_rev_int(k_j, sigma);
    x = fsolve(k_cool,5, options_BE)
    kjInd_FR_rev(s,1)=x;
    save('kjInd_FR_rev.mat', 'kjInd_FR_rev');
end

save('kjInd_FR_rev.mat', 'kjInd_FR_rev');

kappa1 = zeros(nsect,1);
kappa2 = zeros(nsect,1);
theta = zeros(nsect,1);

for s = 1:nsect
fun3 = @(z, k) (1-(z.^2)).*(z.*(exp(z))).^k .*exp(z);
Integral3 = integral(@(z)fun3(z, kjInd_FR_rev(s,1)),0,1);

kappa1(s,1)=kjInd_FR_rev(s,1).*exp(-kjInd_FR_rev(s,1)-1).*Integral3;
kappa2(s,1)=kappa1(s,1)./kjInd_FR_rev(s,1);
theta(s,1)=(kappa2(s,1).*(kjInd_FR_rev(s,1)+1))./((1./(kjInd_FR_rev(s,1)+1))-(kappa1(s,1)+kappa2(s,1))  );
end
thetajInd_FR_rev=theta;
save('thetajInd_FR_rev.mat', 'thetajInd_FR_rev');


init = ones(nsect,1);
CC=Data_FR_revenue(:,4);
AA=theta;
ns=s;
beta_cool= @(B)eqsystem(AA, B, CC, ns);
    x = fsolve(beta_cool, init, options_BE)
    betajInd_FR_rev=x;
    save('betajInd_FR_rev.mat', 'betajInd_FR_rev');
    
